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Abstract 

A new algorithm for reconstructing a two dimensional object from a set of one 
dimensional projected views is presented that is both computationally exact 
and experimentally practical. The algorithm has a computational complexity 
of 0(nlog2 n) with n = for an N x N image, is robust in the presence of 
noise and produces no artefacts in the reconstruction process, as is the case with 
conventional tomographic methods. The reconstruction process is approximation 
free because the object is assumed to be discrete and utiHses fully discrete Radon 
transforms. Noise in the projection data can be suppressed further by introducing 
redundancy in the reconstruction. The number of projections required for exact 
reconstruction and the response to noise can be controlled without comprising the 
digital nature of the algorithm. The digital projections are those of the Mojette 
Transform, a form of discrete linogram. A simple analytical mapping is developed 
that compacts these projections exactly into symmetric periodic slices within the 
Discrete Fourier Transform. A new digital angle set is constructed that allows 
the periodic slices to completely tile all of the objects Discrete Fourier space. 
Techniques are proposed to acquire these digital projections experimentally to 
enable fast and robust two dimensional reconstructions. 

Keywords: Discrete Radon Transform, Mojette Transform, Discrete 
Tomography, Image Reconstruction, Projection Mapping, Discrete Fourier Slice 
Theorem 



1. Introduction 

A common problem in science is to determine the internal structure of an 
object from indirect measurements. Tomography is concerned with the recovery 
or reconstruction of the internal striicturc of an object from its projected "views" 
or projections. In practical terms, projections provide density profiles of an 
object / at different angles 9, either from absorption or transmission of, for 
example, x-rays. The profiles are directly proportional to the attenuation of 
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the x-ray intensity within the object. This attenuation is proportional to the 
magnitude of the line integral 

Mt)= I f{x,v)M, (1) 

J Le t 

at each translate t and where x^y are two dimensional (2D) spatial coordinates. 
The curves Lg t (usually) form a set of parallel lines t S M (with compact support 
D) within the projection ^e{t) having the measure M on [1]. 

An exact reconstruction of the object can then be made when the set of 
projections fj. is given as 

= :0e [0,^)}, (2) 

where the functional notation for ^g{t) has been dropped for convenience and 
it is assumed that — Ms+tt for simplicity^. Equations (1) and (2) define 
the 2D Radon Transform (RT) of the object first developed for ray lines by 
Radon [2] and Funk [3] when using "great circles". Note that the set (2) is 
uncountably infinite, therefore requiring an infinite number of projections for an 
exact or unambiguous reconstruction of the continuous function / [4]. Thus, the 
reconstruction is always ill-posed when using the RT, since only a finite number 
of measurements can be made practically [5] . 

The conventional solution is to acquire a large number of projections (usually) 
at evenly spaced angles. This has practical implications on the quality of the 
reconstructed image [6]. The uniqueness of the reconstruction directly relates to 
the artefacts known as Ghosts or "invisible distributions" [7]. Ghosts are always 
present in reconstructions when using the methods based on this classical RT [8] . 

Discrete RTs attempt to resolve this issue by assuming that the object 
is partitioned into discrete elements, so that an exact reconstruction may be 
computed in the discrete domain with a finite number of projections. The inverse 
problem of RT is then no longer ill-posed and therefore free from reconstruction 
artefacts. Notable discrete RTs include the arbitrary curve block circulant 
discrete RT [9, 10], the d-Hnes discrete RT [11, 12] and the Fast Slant Stack [13]. 
A short review of these methods, which are not required for this work, is 
provided in AppendixA. However, most discrete methods tend to suffer from 
either experimental limitations or high computational complexity. This paper 
presents a discrete approach that is both experimentally and computationally 
viable that utiHses the Mojette Transform (MT), a discrete RT first constructed 
by Guedon et al. [14]. 

1.1. Mojette Transform 

Let the continuous object / be approximated digitally within a point-lattice 
or digital array A (see Fig. 1). This is a reasonable approximation since the 
reconstruction will eventually be represented and processed digitally. Also, let 
the rays form part of a parallel beam, in order to simplify what follows. Then, 
analogous to Eqs (1) and (2), the projections within the MT are defined as 

^^o„{t)= J2 /(^-y)' (3) 



discussion of the effects of scattering and beam hardening in polycfiromatic x-ray beams 
is beyond tfic scope of this paper. 
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Figure 1: An example of a continuous object (I) and its digital representation (II). The digital 
version is deliberately represented in low resolution for illustrative purposes. Note that the 
object and its underlying array need not be rectangular. 

where dpq ~ tan^^(9/p) (or the vectors [(7,p]), with respect to the array A so 
that p,q £2, and the gcd(p, 9) = 1 (i.e. p and q are coprinie). The fractions t/p 
are known as Farey fractions and are also related to the distribution of prime 
numbers [15, 16]. These rational angle projections have a form equivalent to 
that of the discrete linogram, i.e. the number of bins B is dependent on the 
angle 9pq as 



for each projection of a rectangular P x Q image. This has implications on the 
choice of detector resolution and geometry, which are also discussed later in this 
work. 

The definition of the lines fall into one of two classes, depending on the 

desired sampling being fully discrete or as intensities distributed over the area 
between lattice points. In the Dirac Pixel Model (DPM), the pixel is summed 
to its corresponding bin if and only if the line passes through the centre of the 
pixel. The lines Ft^g then form a set of non-periodic and parallel discrete lines 



with f G Z of an object having convex support. Convex support simply means 
that space does not have any local singularities. In interpolated or fractional 
MTs, the pixels are summed according to the length of the pixel intersected 
by the lines Ft.e^^ and can be used to approximate continuous objects [17]. 
Arbitrarily complicated pixel models are also possible using splines of high order. 
The DPM is assumed for the duration of this paper. Fig. 2 shows a simple 
example of a MT for a 4 x 4 image using three projections. 

The MT may be inverted or reconstructed exactly according to one of two, 
albeit practically unsatisfactory, definitions of the Mojette projection set. Firstly, 
one acquires the projections at all possible angles Opq for a given image size, i.e. 
Opq is the set of all rational vectors of order N in all octants of the half-plane for 
an X image. This definition of the MT can be inverted exactly by smearing 
(or back-projecting) all the projections within this Mojette set [19]. Although 
the number of projections are finite in number, in contrast to the RT, it still 



i?= |p|(Q-l) + |<7|(P-l) + l, 



(4) 




(5) 
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Figure 2: An example of a Mojette Transform for a discrete image of size 4x4 using the 
three projections [1,1], [1,-1] and [1,-2]. The bold lines within the right-hand grid shows a 
possible reconstruction path using a corner-based inversion method [18]. 

utilises a large number projections since O^q is large for modest N . For example, 
the inversion of a 128 x 128 image requires a total of 20088 projections and is 
not computationally viable. 

An alternative and more flexible definition of the Mojette projection set was 
proposed by Katz [6]. Since one needs only as many equations (or translates in 
this case) as pixels in order to reconstruct exactly, Katz pointed out that a set 
satisfying 

M-l M-\ 

P^Y.\V^\ (6) 

where M is the number of projections of a P x Q image, is theoretically sufficient 
to obtain an exact reconstruction. Normand tt al. [20] extended the concept to 
discrete objects within arbitrary convex regions. Although this new definition 
substantially reduces the number of projections required for an exact reconstruc- 
tion, it requires new and more sophisticated inversion methods, as the original 
back-projection approach is no longer applicable. A number of schemes have 
been proposed, including a Conjugate Gradient method [21] and a Geometric 
Graph approach [18]. The former is robust in the presence of noise but is not 
suitably convergent, i.e. an appropriate pre-conditioner is yet to be found, and 
the latter is very sensitive to noise. 

The MT was initially constructed to aid work that modelled the human 
vision system [14]. Since then, applications of the MT include image encoding 
and network transmission [20, 17], as well as work on Asynchronous Transfer 
Mode (ATM), data integrity [22, 23], packet networks via an n-dimensional 
MT [24], lossless networking [25] and in scalable multimedia distribution [26]. 
Image watermarking using a Fourier-based method and the MT has also been 
developed [27, 28] as well as another scheme for crypto-marking and transmission 
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over the Internet of medical images [29, 30]. Finally, the MT has been appHed 
to Computerised Tomography (CT) [31, 19, 32, 33]. See Guedon et at [34] for a 
thorough review of the MT. 

This paper presents a new fast digital algorithm for the MT that is both of 
low computational complexity and robust to noise present in projection data. 
The algorithm, which is denoted as the Fast Mojette Transform (FMT) for 
convenience, symmetrises the MT and its projection set to simplify its inversion. 
A novel analytical and exact mapping of Mojette projections to the Discrete 
Fourier Transform (DFT) is also constructed. This mapping, together with a 
new Mojette set definition, allows the direct and exact "reformatting" of Mojette 
projections so as to completely tile all of its corresponding Discrete Fourier space. 
The set, which ranges from + 1 to N + number of projections for an 
N X N image, completely fills DFT space without interpolation and filtering, thus 
allowing an exact reconstruction. If any filtering is required, it is exact and does 
not introduce any artefacts or approximations in the reconstruction. The issues 
associated with detector geometry and resolution of the Mojette projections (as 
given by Eq. (4)) are also addressed. The result is that the projection angles are 
not evenly spaced, the acquisition of projections takes advantage of the digital 
nature of the reconstruction, the inversion is fast and there are no reconstruction 
artefacts. Section 2.4 proposes an experimental set-up for utilising this FMT, 
while preserving its simple inversion process. The FMT is made possible by 
utilising another discrete RT known as the Discrete Radon Transform (DRT). 

1.2. Discrete Radon Transform 

The Discrete Radon Transform (DRT), which was the first discrete RT 
constructed [35-43], is the RT constructed within the same finite geometry as 
the DFT, i.e. the image is assumed to possess periodic boundaries. The DRT 
provides an exact partition of DFT space in the form of periodic slices, where 
a periodic slice is the one dimensional (ID) DFT of a DRT projection. The 
partitioning is known as the Discrete Fourier Slice Theorem (FST) [39, 43]. The 
fast version of the DRT, that utilises the Fast Fourier Transform (FFT) [44, 45], 
is denoted as the Fast Radon Transform (FRT) in what follows. 

The 2D Cartesian DFT is tiled using ID slices (see Fig. 4) along the congru- 
ences (periodic lines) 

t = y — mx (modiV), (7) 
t = X — psy (modA^), (8) 

for an N X N FFT space with A^ = p" and p is a prime number (hence, including 
powers of two) [46]. The lines (7) and (8) are utilised for the slopes 

m = {m : m < N, m e No} , (9) 
s = {s:s<N/p^se No} , (10) 

and the set of translates t as 

t^{t:t<N, t e No}, (11) 

so that the projections are effectively acquired along the vectors [1, to] and [ps, 1] 
modulo A, while the slices are placed along the vectors [—to, 1] and [1, —ps] in 
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Figure 3: The Discrete Radon Transform (DRT). (a) shows the exact or information preserving 
nature of the DRT. (b) shows the projections within the DRT (dotted lines) and how they are 
numerically mapped to linogram (Mojette Transform) projections Mi (solid lines) using the 
nearest neighbour distance d between pixels sampled. 

FFT space (see Fig. 4). Note that the FRT projections/sHces are all of the same 
size or length, unlike the MT, which makes the FRT symmetric in comparison 
and easy to invert (see Fig. 3(a)). 

To completely tile all elements of the entire space at least once, i.e. cover 
all possible coefficients in the N x N FFT space, a total oi N + ^/p slices (and 
hence projections) are required [46]. For the simplest prime case n = 1 so that 
+ 1 projections are required, the gcd(m, p) = 1 always, so the slices of the 
same translate t may only intersect once at the point (0,t). The space is then 
tiled exactly and exact inversion is possible [43]. For the case n > 1 so that 
N -\- ^/p projections are required, those lines having slopes coprime to N will 
intersect only once at the DC coefficient for a fixed t. However, the slopes not 
coprime to N will intersect p times, leading to under-sampling. The extra ^/p 
projections perpendicular to m, i.e. at slopes s, are needed in order to sample 
all of the image. This results in a certain amount of oversampling, which is 
easily and exactly corrected by dividing each coefficient u for all N + ^/p slices 
by gcd(M, A^) [46]. 

Thus, for a 5f2 x 5f2 object, one needs a total of 768 slices for an exact 
reconstruction via the power of two FRT. The same image would require 522 
slices when padded to the nearest prime number 521 using a fast prime-sized 
FFT via Rader's algorithm [45], which is approximately two to three times slower 
than the power of two FFT (see Fig. 4 for an example of a prime-sized tiling). 
The reconstruction is then simply obtained exactly by using the inverse FFT. 
Reconstruction can be started as soon as slices become available, similar to 
conventional tomographic methods. 

Chandra [47] also showed that the projections can be mapped to the Number 
Theoretic Transform (NTT), an integer-only transform analogous to the FFT 
but with better performance for digital projections. The main issue with the 
FRT is that the projections are taken along the periodic lines (7) and (8), such as 
the slices shown in Fig. 4. This makes the experimental geometry not physically 
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(e) (f) 

Figure 4: The exact but periodic slices within the Fast Fourier Transform (FFT) (or equivalently 
the Fast Radon Transform) for the prime case p = 5. Each colour represents a slice of a 
different slope with the DC coefficient/origin centred (black). Note that each vector shown is 
computed modulo p. (a)-(e) shows the slices with slopes ^5 m < 4 (mod 5) in FFT space, (f) 
shows the row sum (perpendicular) slice in FFT space. The set of p + 1 slices tile all of the 
space exactly once. 
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I N = kN' 



N = kN' 



Figure 5: The redundancy within the Fast Mojette Transform as an N X N image. The 
redundancy parameter k is chosen so that N = kN' is a power of two for a N' X A'^' image. 
The original image is assumed to be positioned at the top left for convenience. 



meaningful, limiting its practical application to tomography. Despite this, the 
FRT has been successfully applied to tomography, but no algorithm has resulted 
that is suitably experimentally and computationally practical [48, 49, 32, 33]. 
The next section describes the first major result of the paper, the fusion of the 
MT and the FRT to produce a discrete RT that is both experimentally and 
computationally viable for tomography. 

2. Fast Mojette Transform 

In this section, the algorithm for fast and practical digital tomography is 
described. Given an N' x N' image, the objective is to exactly tile an x 
FFT space, where N = kN' and A; ^ 1, so as to efficiently reconstruct the image 
with the inverse FFT (iFFT). For the FFT space where fc > 1, redundancy 
provides the mechanism for suppressing the effects of noise within the projections 
(see Fig. 5). 

Assuming that the size of the image is = 2", in order to utilise the 
Cooley-Tukey [44] algorithm for the FFT, let the total number of projections 
M = N + N/2. If each Mojette projection maps to a unique FRT projection or 
slice, then they will tile all of FFT space. Once the mapping is determined to 
be one-to-one, the Mojette projections must be converted to FRT projections. 
A schematic of these processes are given in Fig. 6. It will be shown that both 
these key aspects are possible analytically and are exact in nature. A method is 
also given to generate optimal projection angle sets. It is then shown how one 
could experimentally acquire these projections and apply methods to control the 
noise response within the reconstructions. 

2.1. Finite Projection Mapping 

The fundamental operation for taking Mojette projections is the same as 
the FRT because whole pixel intensities are summed along digital rays and 
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Figure 6: A schematic of the main processes within the Fast Mojette Transform for N = 8, 
The Mojette projections (I), acquired using a new Mojette set definition, are converted to Fast 
Radon Transform projections (II) via a new discrete projection mapping, and placed into Fast 
Fourier Transform space (III) ready for inversion. The structure within the dashed box (II) is 
a logical stop, i.e. the projections arc converted in-place. 



placed into discrete bins, i.e. the DPM. The main difference is periodicity, 
i.e. the MT consist of line segments forming parts of the lines formed from 
the congruences (7) and (8) of the FRT. The segments can be linked using 
the shortest distance d vectors between pixels sampled by these segments (see 
Fig. 3(b)). However, these vectors need to be computed using numerical methods 
which can be computationally inefficient [50]. 

The analytical mapping is as follows. For the case n — 1 in N — p"', the 
projection fj,g^^ of the MT maps to the m projection of the FRT as 

m = pq^^ (modiV), (12) 

where is the multiplicative inverse of q that can be computed easily via the 
Extended Euclidean algorithm. For the case n > 1, the projection fig^^ maps to 
the m and s projections as Eq. (12) when the gcd(p, A^) > 1, and 

2s = p-^q (modiV), (13) 

when the gcd{q,N) > 1, respectively. These mappings cover all possible vec- 
tors/angles because the gcd(p, g) = 1 by definition of the angles. The mapping 
results from solving mq = p (mod N) and 2sp = q (mod N) respectively. Ex- 
amples of how each of these congruences are determined is shown in Fig. 7. 

It is crucial to ensure that the Mojette set has a one-to-one correspondence to 
the FRT set, so that the FRT (and equivalently FFT) space is filled completely. 
This new Mojette set definition, which is finite in number and spans the all 
octants in the half plane, while providing a discrete coverage the interval [0,7r), 
is constructed in the next section. 

2.2. Angle Set 

In this section, an FMT set is constructed that optimally fills FRT space 
completely and avoids Mojette projections that are "long", i.e. projections with 
a large number of translates B, by minimising Eq. (4). This is equivalent to 
minimising the ^i-norm of the Farey fractions/ vectors as 

h = bl + 111 (14) 
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Figure 7: A visual presentation of the mapping of Mojctte Transform projections to tlie Fast 
Radon Transform. Solving for m in (a) and s in (b) results in the mappings (12) and (13). 



for each projection in the Mojette set. The simplest way to achieve this is to 
generate angles through the vectors [i,w] with w £ Z in all octants of the half 
plane. The m or s value is then checked with the FRT set in question for each 
vector [l,w] using the mappings (12) and (13). An example of this set for an 
8x8 image is given Fig. 8(a). Note that this set is not the same as the FRT 
vectors [1, m] and [ps, 1], but shares their simplicity. This angle set is a subset of 
the possible minimal ii vectors and is computationally inexpensive to generate. 
Hence, the set is well suited when adapting the detector or image geometry. 
To generate the minimal £i angle set, one computes 



93 



91 +Pi 



N 



P2 



92 - 91, 



P3 



qi + Pi + N 

P2 



P2 -Pl, 



(15) 



where [ • J is the floor (round-down) operators, beginning the computation with 
[91, Pi] = [0,0] and [92,^2] = [1,-/V] until [93, Ps] = [1,1]- Eq. (15) is a modified 
form of the standard way to generate Farey fractions. An example of this set 
for an 8 X 8 image is given in Fig. 8(b). Once generated, the set must be sorted 
by ascending £1 and the first N + ^/p mappings to the FRT angle set chosen. 
However, the set requires the generation of a large number of fractions and 
sorting of these fractions, making the computation expensive. Thus, this set is 
optimal for an unchanging detector and image geometry, since it need only be 
computed once. 



2.3. Finite Conversion 

Once the projection set is known for a given image and FFT space, the 
Mojette projections need to be converted to the finite FRT projections. This can 
be done by equating Eqs in (5), as well as Eqs (7) and (8), with the mappings 
(12) and (13). For a Mojette translate tjn and an FRT translate t-ji, Eq. (5) 
becomes 

9"^ tM = q^^qy-p q^^x 
= y-pq^^x, 
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(a) (b) 

Figure 8: (a) shows an example of a simple angle set, formed from the vectors [l,a] with 
a = 0, . . . , Af — 1 in all octants in the half plane, used within the Fast Mojette Transform 
for 8x8 array. Note that this set is not necessarily symmetric in all octants, (b) shows an 
example of a true £i minimised angle set for the same array. Both mappings tile or cover FRT 
(and hence FFT) space for the image exactly. The origin is marked as a red point. 



when multiplying by g ^. This expression is equivalent to Eq. (7) via Eq. (12) 
when taken modulo N. Likewise, for the s finite projections 



p ^ tM = p ^qy-pp ^ 



X 



p ^qy 



which is equivalent to Eq. (8) via Eq. (13) when taken modulo N. Similar 
equations result for the angles in the other octants. Hence, 

^ ^{q-^ tM i^oAN), if gcd(p,7V)>l 
^ \p^^ tM (modiV), if gcd{q,N)>l' ^ ' 

The conversion can be done prior to the inversion of the Mojette projections, after 
having been acquired and stored, or as the projections are being acquired. The 
latter is more desirable as the Mojette projections are effectively "compacted", 
since the FRT projections are generally shorter than the Mojette projections 
and all have the same size. The next section discusses experimental methods to 
acquire Mojette projections. 

2.4- Mojette Acquisition 

Computationally, the Mojette projections can be acquired trivially using the 
Eqs in (5), since each pixel (x, y) can be summed to its respective translate t [14]. 
The major result of this paper is the experimental data acquisition geometry, 
which is the topic of the remainder of this section. 

The MT can be acquired experimentally in one of three ways. The simplest 
method is to use a ID detector array of size Bmax, which is typically at most N'^ 
in size for rectangular array A when using a minimal £i angle set (see Fig. 9(a)). 
This detector has to have the maximum resolution r of 

J-^ , (17) 

VP + q 
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(b) 

Figure 9: Illustrations of the proposed experimental set-ups for the Fast Mojette Transform. 
Experiment (a) utilises a large detector with a resolution equal to that of the Mojette projection 
requiring the highest resolution. Experiment (b) utilises a small fixed resolution detector that 
is scanned along at an angle (see inset and Eq. (18)) with respect to the view angle of the 
Mojette projection. The angle (f> is used to ensure that the (fixed width) pixels of the detector 
encompass the Mojette ray sums. 
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for unit pixel area, which is effectively the separation between the translates t 
of the angle 9pq corresponding to i3max (see [51, pg. 20-21] for proofs). Since 
the remaining projections require a lower detector resolution, these projections 
become a sub-region within this detector and can be recovered by sub-sampling 
measurements on the detector. The other redundant parts of the detector in 
these cases may be used to estimate noise on those projections. The projections 
can then be converted exactly to FRT projections using the Eqs in (16) and 
reconstruction process may be started immediately as slices become available. 

For detectors smaller than N"^ in size, the detector can be translated or 
scanned across the object in the direction perpendicular to the projection angle, 
similar to the linogram approach [52] (see Fig. 9(b)). To overcome limitations in 
fixed detector resolution and fixed pixel width or size r, the detector is scanned 
at angle ^2 + (p with respect to 6pq. The angle (p is given as 



for each vector [q,p], where ^/y/p^+q^ is the separation between the translates 
t (see also inset of Fig. 9(b)). This ensures that the entire width of the pixel 
in the detector encompasses its corresponding ray sum without the need for a 
multi-resolution detector. However, the detector sensitivity may be adversely 
affected by the angle (j). 

The final method is to use a multi-resolution detector, either of size Bmax 
(being A'^^ in size) or smaller to scan the detector at the angle ^/2 -I- 9pq in relation 
to the object. When utilising the DPM for a continuous object, the projection 
needs to be convolved with a beam-spread function to account for fractional 
sampling of the array A [53]. For the simple angle set, this function is always 
a triangle with a base of 2w and a height of one for the vector [l,w]. In the 
next section, the issues regarding the suppression and the effects of noise are 
discussed. 

2.5. Noise 

The response of the FMT algorithm to inconsistencies in projection data, 
such as noise from detectors, can be minimised in two non-exclusive ways. Firstly, 
the redundancy parameter k can be used to "spread out" these inconsistencies 
over a larger image area. For noise that is statistical in nature, and therefore 
the same on average for different image sizes, the same noise is present over a 
larger area for larger image sizes. Thus, this noise is not as prominent over the 
reconstruction sub-region. Equivalently, the noise is averaged out due to having 
more translates in the projection data. However, since the larger image area 
requires more projections, it exposes the object to potentially more radiation. 
This is a non-issue for data encoding and/or transmission applications. 

Secondly, the redundancy in the entire image, after the inversion of the 
MT, can be used as part of a iterative algorithm for reducing the noise. For 
example, an ii or ^2-iiorm minimising (iterative) algorithm (applied to the pixel 
values) can be used to make image values in the redundant region converge to 
zero [54, 55]. The result of this minimising algorithm gives an approximation 
to the noise in the projection data, reducing the effects of these inconsistencies 
within the image. This approach will be most useful in very noisy or incomplete 
projection data. 




(18) 
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Figure 10: A noise response numerical simulation of the inverse Fast Mojette Transform 
using the simple angle set. (a) shows the original 128 X 128 image of Lena, (b) shows the 
reconstruction with k = 2 and 3% Gaussian noise, (c) shows the actual errors between (a) and 
(b) as a surface plot. On average, the error is approximately 5 out of 256 grey-scales per pixel. 



The redundancy within the image area can also be utilised to reduce the 
number of projections as described by Chandra et al. [32] for the noise-free 
projection data. Any missing FRT projections cause finite ghosts superimposed 
on the image that can be removed exactly via a non-iterative algorithm applied 
to these ghosts in the redundant area. The number of projections is then just the 
number of rows in the initial object, which is generally much less than N -f ^/p. 
Ghost removal in the presence of noise is still an active area of research. The 
next section describes the performance and noise response of the FMT and the 
MT. 

2.6. Results 

Numerical simulations of the FMT were conducted using the DPM of a 
128 X 128 image of Lena on a Intel™ Core2Duo™ E6600 (2.4GHz) processor. 
Fig. 10 shows an example of a simulated reconstruction using the FMT of the 
image with 3% Gaussian noise present in the projections within a 256 x 256 
FFT space (i.e. the redundancy parameter of k = 2). The reconstruction 
(Fig. 10(b)) shows that the result is stable to moderate levels of noise, having a 
Peak-Signal-to- Noise (PSNR) of approximately 35dB, possibly suitable for lossy 
image and video encoding. The errors present on the reconstruction are shown in 
Fig. 10(c), which has a Root Mean Squared Error (RMSE) of 4.8 grey-scales out 
of a possible 8-bit (256) grey-scales. Note that the noise is uniformly distributed 
and no image related artefacts are apparent. 

Graph 11(a) shows the noise response of the two different angle sets as a 
function of the redundancy parameter k. The parameter values k are defined 
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Figure 11: A comparison of the angle sets for the Fast Mojette Transform (FMT) with a 
128 X 128 image of Lena, (a) shows the Root Mean Squared Error (RMSE) with increasing 
redundancy k. (b) shows the computation times of the inverse FMT and the forward Mojette 
Transform (MT) on a log scale with the time in micro-seconds {/i sees). 
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as powers of two to ensure that the FFT space is also a power of two for the 
Cooley-Tukey algorithm [44]. Increasing the parameter k spreads out the noise 
as predicted and exponentially decreases the error on the reconstruction. The 
simple angle set, which is a subset of the £i minimum set, has a better noise 
response because it possesses more translates than the full ii minimum set. 

Graph 1 1 (b) shows the computation times of the MT and the inverse FMT 
in micro- seconds. The inverse FMT, although an order of magnitude slower than 
the MT for very large images, is suitably fast as the time for reconstructing a 
4096 X 4096 image from 6144 projections is approximately 13 seconds. The time 
could be greatly reduced using Graphical Processing Units (GPUs). Considering 
both graphs, the optimal value for the parameter k appears to be around k = 2 
or 4, in order to balance noise response with the speed of computation. Further 
work on the FMT needs to be done in selecting optimal experimental geometries 
and conducting tomographic experiments with real data. 

Conclusion 

A noise tolerant algorithm for fast digital-to-digital tomography was con- 
structed (see graphs in Fig. 11). The angle set for the Mojette Transform was 
redefined to one that was more symmetric and minimal in size, so that inversion 
could be computed with low complexity (see Eq. (14)). Two choices for the angle 
sets were constructed that were either computationally or experimentally suitable 
for variable or fixed image geometry respectively (see Eqs (12) and (13)). A new 
analytical mapping of these projections was constructed that allowed them to be 
compacted directly and exactly into a Discrete Fourier space of desired size (see 
Eq. 16). Once the projections are mapped, the (robust to noise) reconstruction 
can be obtained with a computational complexity 0{n log2 n) with n — iV^ for 
an iV X TV space (see Fig. 10). Redundancy within this space was used to control 
inconsistencies in projection data, such as detector noise. The redundancy also 
allows the exploitation of £i or ^2-norm minimising (iterative) algorithms to 
values within the space to further reduce noise. Experiments and techniques to 
apply this Fast Mojette Transform for real measurements were also proposed 
(see Fig. 9). 
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AppendixA. Discrete Tomography Review 

Discrete reconstruction methods take advantage of digital image geometry of 
the reconstructions to allow exact inversion schemes. The methods of Beylkin [9] 
and, Kelley and Madisetti [10] utilise projections along arbitrary curves and 
reconstruction via Block-Circulant matrices. However, the method requires the 
storage of a large number of matrices and has an unfavourable computational 
complexity of 0{n^) for a 2D image. 
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Figure A. 12; The d-lines for a P X Q image used to approximate continuous lines [11, 12], 
A line, parametrised by x-intercept h and rise s, is divided into two lines comprising the 
to P/2 — 1 and P/2 to P — 1 pixels respectively. The final y-intercept of the first becomes 
/i + [|J and the x-intercept for the second is redefined as h+ [^^J- 

Gotz and Druckmiiller [11], and Brady [12] independently introduced the 
concept of d-lines, which approximate hnes of arbitrary (including irrational) 
slopes (see Fig. A. 12). Sums along d-lines approximate the classical linogram 
and an iterative algorithm has been developed which can recover the image with 
desired accuracy by Press [56]. 

The last approach is to utilise the pseudo-polar Fourier Transform (FT) [57]. 
The Fast Slant Stack method bypasses the need for the filtering within the FST 
via this pseudo-polar FT having a linogram mapping of the slices [13]. However, 
the coefficients are still interpolated (via the Chirp-Z transform) to Cartesian 
Fourier space in order to use the FFT (see Fig. A. 13). 
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2D two dimensional 2 

RT Radon Transform 2 

CT Computerised Tomography 5 

FT Fourier Transform 17 

FST Fourier Shce Theorem 5 

DRT Discrete Radon Transform 5 

FRT Fast Radon Transform 5 

MT Mojette Transform 2 

FMT Fast Mojette Transform 5 

FFT Fast Fourier Transform 5 

DFT Discrete Fourier Transform 5 

NTT Number Theoretic Transform 6 

DPM Dirac Pixel Model 3 
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